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EQUATIONS  FOR  THE  SIMULATION  OF  PLANET  SATELLITE  ORFIT  DETERMINATION 


1.0  INTRODUCTION 


The  equations  defined  herein  arc  for  a  simulation  of  the  process  of  planet 
artificial  satellite  orbit  parameter  estimation  based  on  DSIF  doppler  measure¬ 
ments  and/or  measurements  of  the  direction  of  local  vertical  from  the  satellite 
to  the  planet.  The' simulation  va a  designed  for  analysis  of  the  problem  of  the 
determination  of  satel  Ite  orbits  about  planets  other  than  Earth,  but  is  appli¬ 
cable  with  minor  modifications  to  Lunar  satellites  as  well.  The  simulation 
consists  of  two  parts.  First  a  set  of  observations  are  computed  which  are  with¬ 
out  error  in  the  physical  model  assumed.  Then  the  model  parameters  are  disturb¬ 
ed  in  the  form  of  an  initial  estimate  and  the  process  of  improving  the  estimate 
in  a  specified  subset  of  the  parameters  from  a  set  of  disturbed  observations  i3 
simulated.  Tbe  method  of  estimation  consists  of  sequential  differential 
corrections  of  the  estimate  of  the  parameter  subset  with  each  observation,  using 
the  Kalman  formulation  of  least  squares  filtering.  The  total  set  of  parameters 
which  may  be  estimated  are: 

4  * 

.  Satellite  Keplcrian  orbital  elements 

.  Mars  and  Sun  primary  gravitational  constants 
.  Planet  oblateness  constant 

,  Solar  radiation  pressure 

.  Mars  velocity  and  position  components 

.  Earth  velocity  and  position  components 

The  physical  ms-del  is  considerably  simpler  than  that  of  the  real  world  becauoc 
the  objective  was  to  design  an  efficient  tool  for  analyzing  the  effects  of  the 
most  significant  error  sources.  Tbe  principal  simplification  cf  the  model  is 
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the  propagation  of  Earth  and  planet  orbita  as  heliocentric  conics,  starting  with 
an  arbitrary  point  out  of  the  table  of  epheraerides  in  the  neighborhood  of  run 
start  time.  Since  the  program  allows  estimation  of  Earth  and  planet  ephemerides, 
this  is  not  as  much  of  a  simplification  as  it  would  at  first  appear.  For  a 
planet  such  as  Mars  with  no  large  satellites,  this  method  might  well  be  applied 
without  modification  in  an  operational  orbit  determination  program.  For  Earth’s 
orbit,  the  Moon-Earth  system  must  be  propagated  in  an  operational  program.  The 
second  model  simplification  is  the  neglect  of  all  terns  in  the  expansion  of  the 
planet  gravitational  potential  except  those  for  mass  and  oblateness  with  the 
assumption  that  the  planet  axis  of  symmetry  is  known.  The  third  important 
simplification  is  that  the  D6IF  measurements  are  taken  to  be  range  rate  measure¬ 
ments,  thus  by-passing  a  number  of  steps  in  the  process  of  reducing  the  doppler 
data. 

The  orbit  of  the  satellite  relative  to  the  planet  is  propagated  by  the  method 
of  variation  of  Keplerian  parameters.  The  time  rate  of  change  of  each  parameter 
due  to  the  assumed  disturbing  functions  is  integrated  to  obtain  the  elements  of 
the  osculating  ellipse  versus  time. 

An  important  advantage  gained  by  estimating  Keplerian  parameters  instead  of 
state  vector  is  the  simplification  cf  the  problem  of  isolating  the  initial 
ambiguity  in  orbit  plane  orientation  about  the  line-of -sight  from  Earth  in  the 
case  when  the  estimation  of  parameters  is  based  on  DSIF  doppler  data  only.  The 
method  also  has  an  advantage  in  the  estimation  process  in  that  the  matrix  of 
first  order  coefficients  of  error  in  the  Keplerian  elements  at  one  point  in 
time  with  respect  to  error  in  the  elements  at  another  point  in  time  is  the 
identity  matrix. 
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Table  1  lists  the  nomenclature  for  coordinate  systems  and  variables  as  used 
in  the  derivation  of  equations.  To  simplify  the  notation  the  time -dependent 
variables  are  generally  vritten without  an  explicit  reference  to  time.  When 
used  in  this  way,  it  is  understood  that  the  symbol  for  the  variable  denotes 
its  value  at  problem  time  t.  When  it  is  necessary  to  relate  the  variables  at 
two  or  more  problem  times,  the  time  cymbol  is  added  as  a  subscript  or  in  normal  v 
functional  notation,  e.g.,  or  C(t).  Subscripts  other  than  t  denote  the 
element  of  a  vector  or  a  matrix.  The  equations  for  the  satellite,  planet, 
and  Earth  ephermeridc3  and  for  the  observations  apply  to  the  propagation  of 
numbers  both  in  the  model  and  in  the  estimate  of  the  model.  A  distinction 
is  made  between  model  and  model  estimate  only  when  necessary  to  clarify  the 
nature  of  the  estimation  process.  The  estimate  of  a  model  parameter  is  de¬ 
noted  by  a  caret  above  the  symbol,  e.g.,  the  estimate  of  :<2  Is  represented  by 

A 

X2. 

i 

The  simulation  program  is  documented  in  DB- 84082-1. 
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TABLE  1 

DEFINITION  OF  SYSTEM  VARIABLES  AND  SYMBOLS 


VARIABLE  SYMBOL 


.  Earth  Equatorial  Cartesian  Reference  Axes 

X1 

.  Direction  of  Vernal  Equinox,  E1950 

X11 

.  Earth  Axle  of  Rotation 

X13 

Unit  vectors 

•  *i3  x  *11 

*12 

♦  Planet  Equatorial  Cart>;3ian  Reference  Axes 

*2 

.  Intersection  of  Equatorial  Plane  and 

x  1 

21 

and  Plane,  of  the  Reference  Meridian 

.  Axis  of  Symmetry 

x23 

Unit  vectora 

'  x23  X  X21 

*22 

.  Planet-Centered  Plane -of -the -Sky  Cartesian 

Reference  Axes 

*3 

.  Planet -Earth  Line  at  Arbitrary  Epoch 

x33 

.  Interaction  of  Planet  Equatorial  Plane 
and  P-ane  of  the  Meridian  through 

• 

X31 

Unit  vectors 

•  *33  *  *31 

X32 

.  Unit  Vector  In  Direction  of  Canopus 

CX 

Unit  vectors 

.  Unit  Vector  Along  Vehicle  -  Sun  Line 

SX  J 

.  Vehicle-Centered  Cartesian  Reference  Axes 

*5 

.  SX  * 

X51 

, 

.  SX  x  CX 

x52 

Unit  vectors 

•  *52  *  s* 

X53  , 
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TABUS  1  (Cant’d) 

VARIABLE 

SYMBOL 

DIMENSION 

.  Satellite  Position  -  Velocity  Vector 

.  Relative  to  Earth  in  reference  system 

n 

6 

.  Relative  to  Planet  in  x0  reference  system 

X2 

6  | 

.  Relative  to  Planet  in  x^  reference  system 

*3  • 

6 

.  Relative  to  Planet  in  reference  system 

X3 

6 

e.g.,  (X2)t- ,  1,  3  are  the  position 

components  in  x0 

(X2)y,  ( »  4,  6  are  the  velocity 

components  in  x0 

k 

.  Satellite  Keplerien  Orbit  Parameters 

Semi-major  axis 

a 

Eccentricity 

e 

♦Argument  of  Per  laps  if* 

co  or 

1  2 

Longitude  of  Ascending  Node 

Ai  otX S 

Inclination 

l'lCrl'2 

Time  from  Last  Periapsis  Passage 

*p 

.  True  Anomaly,  Radius  Vector  at  time 

v,  r 

.  Alternate  Notation  for  Keplerian  Orbit  Parameters: 

Y1  »  (a,  e,  -A-  t^) 

Y1 

6 

Y2  »  (a,  e,  g,  c^t  t^) 

Y2 

6 

♦Subscript  denotes  cartesian  reference  system: 

1,  Planet  Eq, 

cystem;  2, 

Plane-of-tho-sky  system 
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TABI£  1  (Coat'd) 

VARIABLE 

SYMBOL 

DIMENSION 

.  Time  Rate  of  Change  of  Satellite 
Orbital  Elements,  Y1 

.  Due  to  Planet  Oblateness 

dSm 

6 

♦  Due  to  Sun  Differential  Graf. 

Effect 

DYS 

6 

.  Due  to  Solar  Radiation  Pressure 

DYR 

6 

.  Planet  Heliocentric  Orbit 

.  Position  and  Velocity  Relative 
to  Sun  in  x1  System 

XB1 

6 

,  Keplerian  Orbital  Elements 

YB1 

6 

.  Earth  Heliocentric  Orbit 

.  Position  &  Velocity  to  Sun  in 
x^  System 

XB2 

6 

.  Keplerian  Orbital  Elements 

YB2 

6 

.  Radar  Site  Location  &  Velocity 

*  Geocentric  Inertial  Latitude  and 
Longitude 

As,  Fs 

If  1 

.  Relative  to  Earth  Center  in  x^ 
System 

X5 

6 

.  Relative  to  Mars  Center  in  x_ 

System 

IS 2 

6 

.  Range  from  Radar  Site  to  Vehicle 

R 

Scalar 

.  Range  Rate 

• 

R 

Scalar 
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TABLE  1  (Cont’d) 

•* 

VARIABLE 

SYMBOL 

DIMENSION 

.  Direction  of  Planet  Local  Vertical  in  x_ 
System  5 

.  Azimuth 

AZ 

1 

.  Elevation 

EL 

1 

.  Earth  Rotation  Rate 

WE 

Scalar  j 

.  Earth  Radius 

Scalar 

.  Mar 8  Radius 

Scalar 

.  Gravitational  Constants 

Earth 

UX 

Scalar 

Planet 

U2 

Scalar 

Sun 

U3 

Scalar 

.  Planet  Oblateness  Constant 

J2 

Scalar 

.  Error  Sensitivity  Matrices 

.  Keplcrlan  Elements  Y1  to  State  Vector 
Elements  X2 

SI 

6x6 

.  State  Vector  Elements  X2  to  Keplerian 
Elements  Y1 

S1I 

6x6 

.  Keplerianjlements  Y2  to  State  Vector 
Eleaente  X3 

82 

6x6 

.  State  Vector  Elements  X3  to  Keplerian 
Elements  Y2 

S2I 

6x6 

.  Solution  Parameter  Set  at  Step  a  to 

Solution  Parameter  Set  at  Step  c—  1 

S 

2$  x  22 
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TABI£  1  (Cont'd) 

VARIABLE 

SYMBOL 

DIMENSION 

.  Parameter  Error  Covariance  Matrix 

C 

22  x  22 

.  Geometric  Transformation  Matrices 

•  X1  —  X2 

T3 

6x6 

•  x2  —  x3 

T4 

6x6 

.  x2  — *  x5 

T5 

3  x  3 

.  Derivatives  of  Observations  with  Respect  to 

I  System  Parameters 

Dfu 

.  Observation  Type  Subscript 

.  Range  Rate 

.  Azimuth  of  Planet  Local  Vertical 

.  Elevation  of  Planet  Local  Vert iced 

.  Parameter  Type  Subscript 

i  -  i 

C  >  2 

lm  3 

♦  Satellite  Orbital  Elements 

X2:  (YP)1  -  a 

J  - 

1 

(X2)2  -  e 

J  - 

2 

<«>3 

J  - 

3 

.  -JL  2 

J  - 

4 

J  “ 

5 

-  *p 

J  •« 

6 

.  Mars  Ephemeris  Elements  XB1 

J  « 

7,12 

,  Earth  Ephermerls  Elements  XB2 

J  31 

13,18 

.  Solar  Frdiation  Pressure 

J  * 

19 

.  Planet  Oblateness 

J  “ 

20 

.  Sun  Gravitational  Constant 

J  “ 

21 

.  Planet  Gravitational  Constant 

_ i_=. 

22 
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2.0  COORDDIATE  SYSTEMS 

Four  cartel  lam  coordinate  systems  are  involved  in  the  system  model.  K?»rth 
and  planet  initial  position  and  velocity  relative  to  the  sun  are  given  in 
Earth  equatorial  E.1950  coordinates  (x^).  The  planet  position  and  velocity 
vector  relative  to  Earth  and  the  radar  site  loci  are  computed  in  the  Earth 
equatorial  coordinate  system  and  transformed  to  planet  equatorial  coordinates 
(xg).  The  satellite  orbit  relative  to  Kars  is  propagated  in  the  planet 
equatorial  system.  Radar  range  and  range  rate  are  then  computed  In  this 
system.  The  direction  of  planet  local  vertical  from  the  satellite  is  com¬ 
puted  in  the  3 un-Canopus  reference  system  (x^).  The  satellite  state  and  orbit 
parameters  are  estimated  by  differential  correction  of  the  Keplerian  orbital 
elements  relative  to  a  coordinate  system  called  the  plane-of-the-aky  system 
(x^),  the  equator  of  which  is  normal  to  the  planet-Earth  line  at  epoch.  The 
first  epoch  is  at  run  start  time  and  may  remain  fixed  or  be  advanced  periodi¬ 
cally.  Figures  1-3  show  the  relationships  of  the  four  cartesian  systems. 

The  3x3  transformation  matrix  (TT3)  between  x,j  and  x0  is  fixed  and  is  a 
program  input.  Since  both  xQ  and  are  fixed  inertial  systems,  the  6x6 
position-velocity  transformation  matrix  is  given  by  T3  *  (  ^  The 

transformation  matrices  between  x^,  Xy  and  x^  are  derived  below  in  terms  of 
Earth,  Sun,  planet,  and  satellite  position  vectors. 

2.1  Planet  Equatorial  -  Plane -of- the -Sky  Transform  (T4) 

At  epoch,  let: 

XB1  :  Planet  state  vector  relative  to  Sun  in  x^ 

XB2  :  Earth  state  vector  relative  to  Sun  in  x^ 

X1K  *  XB2  -  XB1  ;  Earlh  state  vector  relative  to  planet  in  x^ 
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FIGURE  X 
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3.0  TRANSFORMATIONS  AMONG  SATELLITE  ORBIT  PARAMETER  SETS 

The  satellite  orbit  is  represented  at  various  points  in  the  r>rogrru&  by  the 
position- velocity  vector  in  x^  and  coordinate  systems,  as  veil  as  by  the 
Keplerian  elements  relative  to  either  Xg  or  x-.  Figures  4  and  5  shov  the  orbit 
geometry  in  the  two  coordinate  systems.  The  equations  for  the  transformations 
from  position-velocity  at  a  specified  time  to  Keplerian  elements  and  vice  versa 
are  given  in  Reference  2  end  will  not  be  discussed  here.  Let  the*r  outputs  be 
represented  by  the  mnemonics  CARCON  and  CONCAR.  For  example, 

Y1  =■  CARCON  (X2,  time) 

X2  *=  CONCAR  (Yl) 

The  output  of  CONCAR  also  includes  true  anomaly,  v;  and  distance  to  planet 
center,  V* . 

The  Keplerian  orbit  parameter  variations  are  propagated  in  the  planet  equatorial 
system  in  vhich  the  pole  is  assumed  to  be  the  axis  of  symmetry.  Since  the  para¬ 
meters  are  estimated  in  the  plane  of  the  sky  system,  it  is  necessary  to  convert 
from  71  to  Y2  before  filtering  the  observations  at  a  given  time,  and  then  to  con¬ 
vert  again  to  T1  to  propagate  the  orbit  to  the  next  observation  time.  This  can 
be  accomplished  with  the  transformations  CARCON  and  CONCAR  as  follows: 


X2  (t;) 

*  CONCAR 

(U  (tj,  w  ) 

A 

A 

X3  (t.) 

*  T4  * 

X2  (tj 

A 

A 

Y2  (tt  ) 

*  CARCON 

(  X3  (t_),  ) 

The  caret  above  the  symbols  denotes  estimated  rather  than  model  parameter  sets. 
After  filtering  the  observations,  a  new  estimate  replaces  the  old  estimate  in  72 
and  this  is  converted  to  a  new  estimate  in  Y1  by  reversing  the  order  of  trans- 
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formation.  This  series  of  transformations  contains  redundant  operations  since 
in  going  from  Yi  to  Y2  and  vice  versa  only  the  angles  which  define  the  orientation 
of  the  orbit  plane  and  the  major  axis  to  the  reference  axes  change.  The  equations 
for  a  direct  conversion  of  the  angles  and  <*  follows.  (NOTE:  The  trans¬ 

formation  from  Keplerian  parameters  to  state  vector  is  still  required  in  the 
observation  equations.) 

Let  A1  and  PI  represent  the  unit  vectors  of1  the  major  axis  and  pole  of  the  orbit 
in  the  x^  coordinate  system  and  let  A2  and  P2  represent  the  same  vectors  in  the 
coordinate  system.  (See  Figures  h  and  5).  The  direction  cosines  of  A1  and 
PI  are: 


(Al)1 

* 

cos  W  ^  * 

COS  Wt  1  - 

sinco^  •  coe  ^  *  sin-A^ 

(A1)2 

B 

COS  u)  l  ' 

cos  Jl ^  + 

sin  •  cos  i  ^  cos^lL 

(ai)3 

- 

sinvii  * 

sin  C  ^ 

(Pl)1 

C 

sln-A^  * 

sin  6  ^ 

(Pl)2 

9 

-cos7L1  * 

sin  i  ^ 

(PI)  3 

m 

cos  t  ^ 

A2 

x= 

TT4  • 

A1 

P2 

a 

Tl’h  • 

PI 

Let  N2  represent  the  unit  vector  in  the  direction  of  the  ascending  nods  in  the 
Xj  system  and  let  N  be  the  vector  which  forms  a  right-handed  system  with  A2 

and  P2.  That  is,  rT  «  A2  x  P2 

(m)1  -  (A2)2  •  (P2)3  -  (A2)3  •  (P2)2 
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(HM)2 

'  *(A2)1 

•  (P2)3  +  (A2)3  ■ 

(P2)1 

(nh)3 

=  (A2)1 

•  (P2)2  -  (A2)g  • 

(P2)x 

We  can  write  four  equations  Involving  the  components 

of  N2  and  the  argument  of 

periapsis,  u^,  as  follows: 

N23  =  0 

P2  *  N2 

”  (P2)! 

•  (S2)1  +  (P?)2 

•  (M2)2  .  0 

N2  •  nn 

“  oinu>2 

-  (t)2)1  •  (NR^  + 

(n£)2  •  (m02 

iy?  •  A2 

c  cos  ^2, 

«  (N2)3  •  (A2)x  + 

(N2)2  •  (A2)2 

Solving  for  (N2) 

j_i  (H2)2,  cosWg  and  sinu we  have 

UJ 

[’  3 in  uJ 

7 

V  2 

i  cos  u)2 

1 

^2 

r<->2  ] 

=  arctan  [wrj 

From  the  definition  of  inclination,  we  have: 

(  2 

«  arc  cos  ((P2)^) 

Let  the  mnemonic 

for  this  transformation  be  CONCON. 

For  example; 

Y2 

CONCON  (Yl,  T4) 

or 

71 

CONCON  (Y2,  T4T) 

It  is  necessary 

later  to  calculate  the  time  rate  of  change  in  Y2  as  a  function 

of  the  time  rate 

of  change  in 

Yl.  For  this  we  may  either  write  explicit 

equations  for  the  partial  derivatives  of  Y2  components  with  respect  to  Y1 

components,  or  we  may  calculate  finite  differences. 

The  later  method  was 

used  as  follows: 

SHEET  26 


\<  J  4  BO.  14  34  £>*•<».  4-4  6 


NUMBER  D2-84073-1 
REV  LTR 


y  re 

315*7; 


(COBCOH  (Yl  +  A(Yl);  )  .  CONCON  (Yl)  )/  4(Yl)„- 


*  -  3,  I*,  5 


<?  Y2 


1,  2,  6 


Then: 


d(Y2) 

dt 


d  (Yl)  .  3  Y? 

dt  o)  71 


We  viU  also  find  it  necessary  later  to  calculate  the  partial  derivatives  of  the 
satellite  position  and  velocity  components  with  respect  to  Keplerian  elements 
and  vice  versa.  Again,  the  method  of  finite  differences  was  used;  for  instance, 

*  (CONCAR  (Yl  4-  A(n)c  )  -  CONCAR  (Yl)  )/  ^(Yl); 


U.O  EARTH  AND  PLANET  EFHEMERIDES 


The  ephemcrides  of  the  planet  and  Earth  are  propagated  as  conic  sections.  The 
initial  conditions  are  in  x1  coordinates  taken  from  a  table  of  ephemerides. 
These  are  converted  to  heliocentric  orbit  parameters  by  CARCON: 

YBT  *  CARCON  (XB1,  TIME  l)  planet 

YB2  *  CARCON  (XB2,  TIME  2)  sun 

Position  and  velocity  at  observation  times  t;  are  obtained  by: 

XBl(t;  )  -  CONCAR  (YB1  4  ) 

XB2  (t4‘)  -  CONCAR  (YB2,  t<;  ) 

When  the  position  and  velocity  components  are  included  in  the  set  of  parameters 
to  be  solved  for  the  estimates  in  the  parameter  3ets  are  denoted  by  placing  a 
caret  above  the  vector  bar.  After  filtering  in  cartesian  coord laates,  the  new 
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J.  M.  Dariby  in  reference  2  develops  a  method  for  obtaining  the  time  derivatives 
of  the  conventional  orbit  elements  in  the  frame  of  planet  synmetry  from  partials 
of  the  perturbing  function  with  respect  to  these  elements.  Hi 3  analysis  is 
applicable  to  any  perturbation  expressible  in  terms  of  radial,  cross-radial 
and  normal  force  components  or  a  disturbing  function  F,  F  being  a  function  of 
the  orbit  elements. 
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Assign  the  following  notation: 

a,  e,W,-A,  i  ,  T  conventional  orbit  elements  in  reference  of  ayanetry 


r,  v 

/* 

F 

T 


n. 

R,B,N 


radius  rector  and  true  anomaly  at  time  t 

planet  gravitational  constant 

disturbing  function 

time  at  peri&psis  passage 

t  »  t  -  T,  where  t  is  present  time 

orbit  period 

mean  motion,  2 IT/ P 

Radial,  crosc-radiol,  and  normal  components  cf  disturbing 
force  at  r  and  xr 


and  define  the  intermediate  variables 

4jJ*  *  +  — /l— 

1  c  -nT  +  *■*>+  vt- 

i  1  *  -n(t-T)  +  *>  +-d- 

U  «  UO  ♦  or 

From  equations  11.9*9  we  select  expressions  for  the  time  derivatives  for  five 
of  cur  aix  orbit  elements. 


da 

dt 

2na^  . 
/+~ 

jt 

de 

na(l-ec 

)  .  ££ 

dt 

/*e 

d(J 

na  T5T 

ef  «JF 

dt 

/e 

<*c 

d  A. 

na 

c)F 

,r~  £S 

dt 

W*-® 

r  c8ecj7 

na  -j  l-e^ 
/ue 


.  <9F 

S%  +  ^ 


cot  <•  ^5 
.  cJc 
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/  o"SL 


ceec  t  +  taxi 


±i{  41  +  il 

2  £  >t.  ^<-0* 


also  includes  an  expression  for  the  time  rate  of  change  of 


-Pna2 .  +  na  Jl-e2 

/*-  <2  a  ,ze 


1  -  y/l-e4 


2  1 


Z<  Jl-e2 


An  expression  for  rate  of  change  of  time  from  pericenter  la  obtained  by 


differentiating  the  equation  for 


!!b  «  -JL.stll  ♦  _= — iisf  .  i._k.«s 

dt  n  dt  n  dc  2  a  dt 


.  i  .  l£i  +  1  .&*  +  I  .if  .  3  .  JE  -d* 

n  dt  n  dt  n  dt  2  a  dt 


Substituting  equation  6  yields  .* 


7)  .  issl'.a  .  ™  r J7_esf^ 

dt  n/*  n^e  1  ^  1  )  Je 


-  - tan  Ai'-i?  +  +  IJbi^  -i.J8.ift 

•$,u,  Ji_e2  2  n  dt  n  dt  2  a  dt 

Danby’s  equations  11-9.8  gives  the  required  derivatives  of  F: 

8>  H  ■  «i 


p  f  r 

^  m  -R  a  oosv  +  Ba  j  j  +1  einV 


T*  *  Nr  sin  u 


p  _ _ 

iz  jgjsj&r  .  B  a  -U-e2  + 
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^  p  2  1  *  * 

3-2)  *  -  By  sin  -g-*  -  Nr  sin  ^  cos  u 

13)  *f  -  —4.  .  sinir  -  £2!  . 

*  h-o2 

Substituting  expressions  8)to  13)  in  expressions  l)  to  6)  yields  the  time  rate 
of  change  of  orbit  elements  due  to  each  disturbing  function  F. 


5.2  Hie  Disturbing  Function  Due  to  Ob lateness 
Hie  oblateness  disturbing  function  is 

F  "  ■  J2  (  •  3-«ln2o  •  aln2  -  1  ) 

where  rQ  is  the  equatorial  radius  of  planet  and  is  the  oblateness  constant. 
The  partials  of  F  with  respect  to  «*  and  vJ  can  be  found  by  straight  forward 
differentiation. 

2 

F  r0  2  /  v  » 

— - —  -  sin  (w  +  v)'sin*  •  cos  t* 


and 


-  -  3/^ Jg  '  *L  *  a*n  ( T")  '  co a  (lO’-f  lt) 


Since  symmetry  about  the  axis  of  rotation  is  assumed  (inclusion  of  only) 

n  -  0 

J-SL 

The  force  components  R,  N,  and  B  can  now  be  defined  In  terms  of  these  derivatives. 
Equation  10)  gives  an  expression  for  N  • 


N 


sj:  .  — 1 - 

r  s*n  1  *>+  w 
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C 


i  o 


Equation  12)  gives  an  expression  for  B  : 


B  -  1  sin  t  cos  ' u  +  v) 

^  ‘  r  sin^  ^ L  sin  ( i-J  +  v') 

Finally,  equation  !l)  gives  an  expression  for  R  • 


R 


since 


-jEh^L  f  jg  .  f .1  ■  i:  g“8.  {  ?.+.  }fl 

•e.invUw  ]t  Sl/X8lr  ^’JV  3) 


-  O 


«*<->  * 


These  values  of  R,  B,  and  N  by  direct  substitution  into  8),  9),  and  13)  give 
isasediate  values  for  the  remaining  parti als  of  F  (oblate)  with  respect  to  the 
parameters. 

5.3  The  Disturbing  Functions  Due  to  Sun  Mass  and  Solar  Radiation  Pressure 

The  position  and  velocity  of  mars  in  sun  centered  earth  equatorial  cartesian 
coordinates  is  denoted  by  XBi.  These  are  transformed  to  sun-centered,  mars* 
equatorial,  cartesian  coordinates,  XS. 


XB 


T3  *  XBI 


Sun 


Satellite 


SX2 

/ 

y 

4- 


\  xs- 


X2 


_ Ni 


Planet 
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The  position  and  velocity  of  the  sun  in  oars  centered  oars  equatorial  cartesian 
coordinates  given  by 

5X2  =  -  X?  -  XS 

The  disturbing  function  due  to  sun  mass  in  satellite  centered  mars  equatorial 
coordinates  is  then 


SF  *  GM 


sun 


XS 


3X2 


v  /xS  |3  fsX2(3 

The  rotation  matrix  from  mars  equatorial  to  radial,  cross-radial,  and  normal 


i6  denoted  by 

tt6. 

TT6(1,1) 

s 

COS  (  +  it)  -  COS  -A 

tt6U,2) 

B 

cos  ( +  ir)-sin_A. 

TT6(l,3) 

S 

6in  ( <*>  *  v")  •  Bln  0 

TTb(?,l) 

= 

-  cos  -A  *sln  (  v~) 

TT6(2,2) 

- 

-  sin  (  u>  +  is)  -  sin  -A 

Tf6(?,3) 

= 

cos  ( v)  •  sin  i 

TT6(3,1) 

B 

sin  l  -sin 

tt6(3,2) 

B 

-sin  i  -cos -A- 

TT6(3,3) 

3 

cos  f 

Then  the  radial,  cross-radial,  and  normal  force  components  are 


V 

SF(l) 

B 

«(tt6)- 

SF(?) 

L  »J 

SF 

_SF(3)„ 

These  can  nov  be  substituted  directly  into  equations  8)  to  13)  to  give  partials 
of  FCsun  mass) with  respect  to  the  orbit  parameters. 
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RF  -  -  RCOH  •  TT6  *  S3( 

where  RCON  is  a  constant  depending  on  vehicle  characteristics.  Again,  substi¬ 
tution  in  equations  3)  to  13)  gives  the  partials  of  the  disturbing  function  with 
respect  to  the  parameters. 

5.4  Extension  of  the  Model 

Tfte  model  for  variation  of  parameters  is  extended  by  defining  each  new  disturbing 

function  and  deriving  its  partied  derivatives  with  respect  to  the  parameter*  in 

one  of  the  two  ways  illustrated  above.  If  the  disturbance  is  represented  by 

one  or  more  terms  in  the  expansion  of  the  planet  gravitational  potential,  the 

parti&ls  are  obtained  by  differentiation  of  the  terms  with  respect  to  the 

% 

parameters  as  was  done  for  the  cblateness  function.  If  the  disturbing  function 
is  expressed  as  a  force  vector  as  in  the  case  of  sun  mass  and  radiation  pressure, 
the  vector  is  resolved  into  components  R,  N,  B,  which  are  substituted  in 
equations  8)  to  13)  to  obtain  the  partials  of  the  disturbing  function  with 
respect  to  the  parameters.  Equations  l)  to  6)  apply  for  all  disturbing  forces 
to  obtain  the  time  ra*e  of  change  in  the  parameters. 


6.0  THE  OBSERVATION  EQUATIONS 


The  range  of  the  satellite  from  a  radar  site  Is  given  by 


where  (X2)  and  (X52)  ,  ( c  *  1,3)  are  the  components  of  the  satellite  and 


u *  4m:  1434  0»>Ci,  4-61! 
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site  positions  relative  to  the  planet  in  the  coordinate  system.  The  site 
position- velocity  vector  is  calculated  first  in  the  system,  then  transformed 
to  aQ  follows: 

r .  -  ( r  »)0  ♦  we  (t  -  to) 


(OS)1 

=  RE  • 

c°a a 

•  cos  r 

l  '  s 

+  (X1E) 

(5QS)g 

=  RE  *  cos  X 

s 

sin  r  + 

s 

(X1E)2 

(X1S)3 

-  RE  • 

3lni  a 

-  +  (xis) 

3 

(jqs)4 

=  -WE  •  RE 

cos^a 

•  sin  r  s 

+  (XUB) 

(-tlS)j 

*  WE  •  RE 

•  cos  ra 

+  (XXE) 

<X“>6 

*  (X.1K)6 

XS2 

T3  * 

xls 

The  filter  equations  require  the  formulation  of  the  partial  derivatives  of  the 
observations  with  respect  to  the  parameters.  In  all  of  thebe  formulations,  to 
be  derived  later,  the  set  of  partied  derivatives  of  the  observations  with 
respect  to  the  elements  of  the  position-velocity  vector  xS  appears  as  a  term. 
Differentiating  the  expression  for  range  rate  with  respect  to  the  elements  of 
X2,  we  have: 
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1*  »»  3 


H 

T S77 


(«)«:.  3  - 


<•  -  U,  t>,  6 


J  D  ,  * 

For  future  reference,  let  ^  *  1,1,  •  **  ,  U 


6.2  Planet  local  Vertical 


The  position  vector  of  the  vehicle  relative  to  the  planet  in  the  Sun-Canopua 
reference  system  is  given  by: 

X5  *  T5  *  >5 

where  T5  is  the  transformation  matrix  derived  in  paragraph  2.2. 


The  direction  of  planet  local  vertioal  is  defined  in  terms  of  two  angles  in 
the  Sun -Canopus  reference  system  as  follows: 


"  (X5)n 

arctan - 

(X5)1 

r  <X5) 

3 

arcsin  » 


Differentiating,  we  have  the  partial  derivatives  of  AZ  and  K1  with  respect  to 
the  satellite  position  vector  X5: 


1U5),  "  (  '(X5)2  •  eo*8{«)/  (X5>1 


^  AZ 

3TxT 


cos  (AZ, 

ihu 


shrei  36 


COMPANY 


NUMBER  Dg.84073-1 
REV  LTR 


■m  - 

(X5)3  • 

(X5)x/  /X 51  • 

cos  EL 

-ast  ■ 

(X5)3  • 

(X5)g  /  1X51  ■ 

cos  EL 

-Mi  ■ 

(1+  «X5) 

/>/'X5f  • 

cos  EL 

Applying  the  chain  rule: 

3>(AZ) 

•  T5 

3(X2) 

^(X5) 

cJ(X^) 

<*(X 5) 

^ (El  ) 

)  . 

S&1  . 

^(AZ) 

•  T5 

^(X2) 

j(X5) 

X2 

•  5  OS) 

For  future  reference,  let 

D1£,J  ‘ 

<5(az) 

3Tx2) . 

m  _ 

Mkl) 

D13,J 

K 

7.0  THE  PARTIAL  DERIVATIVES  OF  OBSERVATIONS  WITH  RESPECT  TO  THE  PARAMETERS 

The  matrix  of  partial  derivatives  of 

the  observations  with  resrect  to  the  para- 

meters  to  be  solved  for  is  denoted  by  D*^,  with  the 

subscript  definitions  as 

outlined  on  the  last  page 

of  Table  1, 

►  The  equations  for  the  partial  derivatives 

of  the  observations  with  respect  to  vehicle  position  and  velocity  elements,  X2, 
were  derived  in  paragraph  6.1  and  6.2.  This  matrix  of  partial©  i3  represented 
by  D1  *  ^  Ca  and  J  *  ly-jC.Dl  and  D  are  related  to  the  matrix  of  partial 
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derivatives  of  X2  elements  with  respect  to  the  parameters  by  the  chain  rule  of 


d If f e  rentlation : 


D  «=  D1 


^ (XP)X 

» 

3Nt 


5(xs)1 


{.<2).  '  ‘  '  (3), 


Ihe  remainder  of  the  problem  reduces  to  deriving  the  equations  for  the 
derivatives  of  X2  with  respect  to  the  complete  parameter  set  Y. 


The  first  group  of  parameters  consists  of  the  Kcplerian  elements  relative  to 
the  plane  of  the  sky,  Y2  •  (Y  •.  =  (Y2);  ,  i  *  1,-6) 


By  the  method  of  finite  differences: 


*  [CONCAR  (Y1  +  ^(Yl)t  )  -  CONCAR(Yl)  ]/  ^(Yl),.* 


=  i,-.6 


Di  »  *  [cONCC'N  (Y2  +  ^YS)^  -  CONCON  (Y2)  J  /  *(Y£)t- 


<•  *  1,-6 


D1  •  D2  •  D3 


<•  »  l»-3 
J  »  1,-6 


The  second  group  of  parameters  consists  of  elements  of  planet  and  Earth  velocity 
and  position  vectors.  The  direction  of  planet  local  vertical  is  insensitive  to 
planet  and  Earth  ephemerls  deviations,  so: 
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Tho  derivative  of  range  rate  with  respect  to  planet  ephemeris  deviations  is 
given  'by: 


^  R 


Xxp) 


X2  ) 

eKxi) 


<Mxi) 

^oa 


or  D 


D1  •  T3 

l  -  1;  J  *  7,«-l? 


‘Hie  derivative  of  range  rate  with  respect  to  Earth  ephemeris  is  given  by: 


d  R 


c)  R 


<X.XB2/ 


Xxss) 


o/_XS2  ) 

<>xis' 


~ ( X1S  1 
(  XB!>J 


where  XS2  and  XIS  are  radar  site  position  and  velocity  vectors  in  Xg  and  x^, 
coordinate  systems,  respectively.  In  matrix  notation: 


•J 


*  -D1 


T3 


1;  J  -  13,-13 


The  final  group  of  parameters  are  the  physical  constants:  the  planet  and  Sun 
primary  gravitational  constants,  planet  oblateness,  and  solar  radiation  pressure 
The  derivative  of  X2  elements  with  respect  to  planet  gravitational  constant  is 
found  by  the  method  of  finite  difference  as  follows: 


^  CGNCAK  (Yl,  Ug  ♦  AUg)  - 


*  U* 


CONCAR  (Yl,  U2) 


A  If- 


Then 


c>  R 


5  R 


or  D(  1,22)  S  D1  » 


Tuk 


^Gcg) 
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The  elements  of  ^  for  the  remaining  physical  constants  are  functions  of 
the  variation  in  the  Keplerian  elements  associated  with  each  constant. 
Consider  the  effect  of  oblateness  over  an  observation  interval  At: 


X  +  At 


Yl(t)  -  Y1  (t  -  At)  + 


dt 

oblate 


Yl  (t  -  At)  ♦  J2 


(AYl)  Oblate 
.72 


t,  t  +  At 


AYl)  Oblato 
J  2 


AYM 

"j r~ 


t,t  ♦  At 


Applying  the  chain  rule  of  differentiation: 


(Q3SKRVATI0No’ 
1  J? 


*  D, 


>(0B3)  .  )(*?)  .  AYM 

2  (X2)  H7I)  J2 


AYM./J2 


1,  2,  3 


Similarly  for  solar  radiation  pressure  and  Sun  mass: 


*=  D1  •  D2  •  [aYR/RCON  AY3/U  1 
L  1 

*  1,  2,  3;  J  *  19  and  21 
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8,0  HIE  KAIXAN  FILTER 

The  general  theory  of  estimation  which  Justifies  the  process  of  least  squares 
filtering  is  treated  in  references  3  through  6,  Particularly,  the  equations  for 
the  Kalman  filter  used  in  this  simulation  are  developed  in  reference  6.  The 
equations  as  implemented  produce  a  differential  correction  in  the  current  estimate 
of  parameters  with  every  observation  set  taken  sequentially  in  time.  The  correct¬ 
ion  is  a  function  of  the  following  variables: 

the  matrix  of  derivatives  at  the  observation  time  of  the  observations 
with  respect  to  the  parameters,  D,  where  t>  is  of  dimension  IXM,  I  being 
the  number  of  observation  types,  and  M  being  the  number  of  parameters  in 
the  set  being  estimated; 

.  the  estimate  of  error  in  the  parameters  before  the  correction  is  made, 
represented  by  the  parameter  error  covariance  matrix  C  of  dimension  Mfll. 

.  the  observation  error  covariance  matrix  J£N  of  dimension  XXX • 

the  vector  of  residuals  in  the  observations  before  the  correction  is 
made,  represented  by  the  vector  F  of  dimension  I. 

. 

The  estimate  of  observations  and  the  associated  elements  of  the  D  matrix  at  time 
t  are  calculated  as  a  function  of  the  estimate  of  parameters  at  time  t  before 

the  correction  is  tdb&e,  using  the  equations  developed  in  paragraphs  2  through  7* 

. 

*ftie  residual  vector  F  at  time  t  +  6  t  is  the  difference  between  the  simulated 

observations  (Model  observations  plus  random  error)  and  predicted  observations 

based  on  the  parau^ter  estimate  at  time  t.  The  matrix  KN  is  a  program  input 

which  is  constant  throughout  the  run.  The  matrix  C  at  time  t,  before  the 

correction  is  made,  is  calculated  as  a  function  of  its  value  after  the  last 

correction  at  time  t  -  d  t,  and  as  a  function  of  the  first  order  deviations  in 

the  parameters  at  time  t  with  respect  to  deviations  in  the  parameters  at  time 

. .  . . . .  . . . . . . . . — . ■  „  — . . 
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t  ~  A  t.  I*t  the  matrix  of  coefficient  deviations  at  time  t  vith  respect  to 
deviations  at  time  t  •  A  t  be  denoted  by  S.  In  the  parameter  set  ve  have 
chosen,  all  of  the  diagonal  elements  of  S  are  one  and  most  of  the  off-diagonal 
elements  are  zero.  The  exceptions  are: 

•  the  coefficients  of  deviations  in  the  Kepler ian  elements  at  time  t 
vith  respect  to  the  ob lateness,  Sun  mass,  and  solar  radiation  pressure 
constants  are  respectively  ZHM/J2,  Ifs/tly  and  25*  JtCQH .  These  are  the 
elements  ^£>y  tm  1  through  6,  and  J  «  19,  20,  and  22. 

•  the  coefficients  of  deviations  in  planet  and  Earth  position  components 
at  time  t  vith  respect  to  velocity  component  deviations  at  time 

t  -  At  are  all  equal  to  dt. 

The  parameter  error  covariance  matrix  at  time  t  before  correction  is  given  by: 

ct  •  3  •  cwt  •  sT 

The  differential  correction  in  the  parameters  at  time  t  is  given  by: 

OTt  -  V  •  Ft 

t  r  v  1  -  j. 

vhere  W  »  Ct  ‘  Dt  I  Dt  '  ct  '  *  £N  J 

The  nev  estimate  in  the  parameters  and  parameter  error  covariance  are: 

KStf  \  -  Yt  +  DYt 

**et  m  Ct  *  W  *  Dt  ‘  Ct 
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REVISIONS 


DESCRIPTION 


DATE  APPROVAL 


Revision  (A)  consists  of  the  following  corrections  and 
additions  to  the  text: 

l)  Equation  (T)  on  page  30  Is  revised  to  drop  the  last 
term  in  the  equation.  Replace  page  30  with  revised 
page  30. 


2-8-66 


m 


The  first  half  of  page  32  Is  revised  to  Include 
equations  vhi?h  eliminate  a  problem  of  discontinuity 
in  the  original  equations  at  true  anomaly  equal  to 
zero.  Replace  page  32  with  the  revised  page  32  and 
a  new  page  32a. 


3)  A  typographical  error  Is  corrected  in v the  last  line 
of  page  36.  Replace  page  36  with  the  revised  page  36. 


h)  Errors  in  the  second,  third,  and  fourth  equations  on 
page  37  are  corrected.  Replace  page  37  with  revised 
page  37. 
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5) 


dt 


-na 


at  "  r-  ~T  { £  ‘  +  *an  ~r  L\H  *  li’  } 

u  / 1  -  c 


He  also  Includes  an  expression  for  the  time  rate  of  change  of 

dIi_  .  *ml  .  m  *  a ,yI7°g  V1  .  /i77! 


;> 


dt 


3a 


p.e 


na 


J  Y 


T  tQn  4" ‘  "If 


An  express  ion  for  rate  of  change  of  time  from  pericenter  is  obtained  by 
differentiating  the  equation  for  1^. 


dt 

— j 

dt 


JL 

•n 


_1 

dt 

_L 

n 


1  .  da>*  _  3  .  *p  .  da 

n  dt  "  2  a  dt 


d^ 

dt 


J__  .  do 
n  dt 


d fl 

dt 


Substituting  equation  6  yields: 


dt 

7)  — -2 
' ;  dt 


2nac 

a* 


d? 

3a 


na 


/T-~?  n — rr  l  dF 

•  ’  —  — ■  4  1  -  / 1  -  e4*  .  — r- 
n  u  e  l  J  de 


na 


n  u  /l  -  e 


X 


tan 


1  •  .  d  F  1  ,  dtr, 

2  ITT  n  dt 


dO 

dt 


Danby's  equations  11. 9*8  gives  the  required  derivatives  of  F: 
„  r 


8) 

9) 


11 

3a 

11 

3e 

d  F 


*  -R  a  cos  v  ♦  b  a 


tt  •>} 


sin  v 


10)  «  N  r  sin  u 


11) 


d  F 
3  uj* 


«P  ae  sin  v  .  B  -A-  •  /{Te5" 

^777“ 


+  Br 


v  *  A 
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SHEEi  30 


NUMBER  E2-84G73-1 

REV  LTR  a 


Substituting  the  expressions  for  dE  from  Equation  (17)  successively  In  Equation 
(1 6)  yields: 


<■«>  #  -  *  • 


sin  v  r 

T~  4  “ 


«  jlL  .  JSUL 
da  P  J 


.  /TZ1F 


where  M  ig  the  aean  anomaly 


Finally, 


at  time  t. 


A1  B  Jl 

o  a  dr 


+  dF 
da 


and  R  ’  Cr)  C*H  ) 


5*3  The  Disturbing  Functions  due  to  Sun  Mass  and  Solar  Radiation  Pressure 

The  position  and  velocity  of  Mars  in  sun  centered  earth  equatorial  cartesian 
coordinate  is  denoted  by  XBI.  These  are  transformed  to  sun- centered,  Mars- 
equatorial,  cartesian  coordinates,  IS. 

XS  *  T3  •  XB1 


\  >l3 


Satellite 


Planet 


V  •  480.'  T4,f4  OHlG,  4-n' 
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Equation  (12)  yeilda  an  expression  for  B: 


B  * 


sin  i  eo3  (u>  *■  w 


x*  2  / 

sin*"  vtv  sin  (iv  +  v) 


The  apparent  discontinuity  at  («u  +  v)  «  0  due  to  the  sin  (jj  +  v)  tens  in  the 
denominators  of  N  and  B  does  not  actually  exist  because  dT/tf  has  a  factor  sin* 
(o>  +  v).  Equation  (ll)  yields  a  solution  for  the  component  R  at  all  points 
except  at  v  *  0.  To  avoid  this  discontinuity, the  partial  of  F  with  respect  to 
the  semi-inajor  axis  is  derived  and  substituted  in  Equation (8)  to  yield  in  ex¬ 
pression  for  R.  Equations  (9)  and  ( 13)  then  yield  the  remaining  partiaia  of  F 
with  respect  to  the  elements. 


The  partiaia  of  F  with  respect  to  r  and  v  are: 


~  t*  J2  ( sin^  c  )  (sin 2  (a  ♦  v)  j 


-  3  *  J 


2  r3 


(sin^v')  ^»in  (ui  +  v)  j  (cos  v)  j 


The  partiaia  of  r  and  v  with  respect  to  (a)  are  derived  from  the  ellipse  tia*  and 
position  equations  as  follows: 


(14)  (t  -  T)  *  ♦  /  —  (E  -  e  sin  E) 


(15)  cos  E 


a  -  r  m  e  +  coa  v 
ae  1  +  e  coo  v 


Holding  time  constant  and  differentiating  (14)  and  (15)  with  respect  to  (a) 
yields: 

(16)  (E  -  e  sin  E)  da  »  -a  (l  -  e  cos  E)  dE 


(IT)  dE 


ae  sin  £ 


•  (dr  -  X  da)  or  dE  ■  i~"H— v  dv 
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For  future  reference,  let  Dlj  ^  «  ■  ( Jgij* .  *  ^  *  1,  2,  ...  6 

6.2  Planet  Local  Vertical 

The  position  vector  of  the  vehicle  relative  to  the  planet  in  the  Sun-Canopus 
reference  system  is  given  by: 

X>  *  T5  *  *2 

where  T5  is  the  transformation  matrix  derived  in  paragraph  2.2. 

The  direction  of  planet  local  vertical  is  defined  in  terms  of  two  angle*  in  the 
Sun-C&nopus  reference  system  a*  follows: 

(X5)?  , 

“  "  arctan  l  THT l  j 

(X5),  , 

EL  -  arc^ln  L  -p^p  J 

Differentiating,  we  have  the  partial  derivations  of  AS  and  3L  with  respect  to 
the  satellite  position  vector  X5: 

M)1  "  V*  (XS)2  ■  C0*2 

d  (A2l  cos2  (AZ) 

ttx5)2  3 
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Applying  the  chain 

rule: 

?(AS) 

ikz 

4C&) 

t  (■<*>) 

.  . 

H<2) 

<)(~ 5) 

'  (g.) 

<■  \xzj 

d  (el) 
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&  >2 
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For  future  reference,  let 

Ei2,J 

d  (AZ) 

"Hi 

C13y 

• 

7.0  THE  PARTIAL  DERIVATIVES  OF  OBSERVATI ONS  WITH  RESPECT  TO  THE  PARAMETERS 


The  matrix  of  partial  derivatives  of  the  observation®  with  respect  to  the  para* 
meters  to  be  3olved  for  is  denoted  by  D,  4,  vlth  the  subscript  definitions  as 
outlined  on  the  last  page  of  Table  1.  The  equations  for  the  partial  derivatives 
,of  the  observations  with  respect  to  vehicle  position  and  velocity  elements,  X2, 
were  derived  In  paragraph  6.1  and  6.2.  This  matrix  of  partlals  is  represented 

by  D1 .  1,2, }  and  j  *  1,..6.  D1  and  D  are  related  to  the  matrix  of  partial 

( to 
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